package math;

import dataTypes.ISequenceData;

/**
 * A helper class with methods to do quality related calculations.
 * 
 * @author Christian Johansen
 * 
 */
public abstract class Quality {
	private static final double PHRED2PCONST = -0.230258509;
	private static final double P2PHREDCONST = -4.342944819;

	/**
	 * Converts a quality value from phred format to probability of error.
	 * 
	 * @param phredValue
	 *            phred value of reading.
	 * @return probability of base reading being an error.
	 */
	public static double phred2p(int phredValue) {
		return Math.pow(Math.E, PHRED2PCONST * phredValue);
	}

	/**
	 * Converts a quality value from probability of error to phred format.
	 * 
	 * @param pError
	 *            probability of base reading being an error.
	 * @return phred value of reading.
	 */
	public static int p2phred(double pError) {
		return Math.round((float) (Math.log(pError) * P2PHREDCONST));
	}

	/**
	 * Calculates the average quality using a number of fixed size sliding
	 * windows.
	 * 
	 * @param seqData
	 *            sequence data.
	 * @param startSize
	 *            smallest window to use (odd number).
	 * @param step
	 *            increment between the window sizes to use (even number).
	 * @param numberOfSizes
	 *            number of different sliding window sizes to calculate.
	 * @return the average quality per base per sliding window size.
	 */
	public static double[][] calcWindowQuality(ISequenceData seqData,
			int startSize, int step, int numberOfSizes) {

		int seqLength = seqData.getSeqLength();

		byte[] seqQuality = seqData.getSeqQuality();

		double[][] qualAvg = new double[numberOfSizes][seqLength];

		int windowSize = startSize;
		for (int i = 0; i < numberOfSizes; i++) {

			for (int seqPos = 0; seqPos < seqLength; seqPos++) {
				int windowSum = 0;
				int currentWindowStart = Math.max(0,
						(int) (seqPos - (int) (windowSize / 2)));

				int currentWindowEnd = Math.min(seqLength - 1,
						(int) (seqPos + (int) (windowSize / 2)));
				for (int pos = currentWindowStart; pos <= currentWindowEnd; pos++) {
					windowSum += seqQuality[pos];
				}
				qualAvg[i][seqPos] = (double) windowSum
						/ (1 + currentWindowEnd - currentWindowStart);
			}
			windowSize += step;
		}
		return qualAvg;
	}

	/**
	 * Calculates the average quality using a dynamic sized sliding window that
	 * adjusts its size so that it covers an area with quality difference >=
	 * minMaxThreshold
	 * 
	 * @param seqData
	 *            sequence data.
	 * @param minMaxThreshold
	 *            the minimum difference between min and max quality in a
	 *            window.
	 * @return the averaged quality.
	 */
	public static double[] calcDynamicWindowQuality(ISequenceData seqData,
			int minMaxThreshold) {
		int seqLength = seqData.getSeqLength();

		byte[] quality = seqData.getSeqQuality();

		double[] toReturn = new double[seqLength];

		for (int i = 0; i < seqLength; i++) {
			int windowSize = 1;
			int max = 0;
			int min = 0;
			while (max - min < minMaxThreshold) {
				windowSize += 2;
				max = max(quality, i - (int) (windowSize / 2), i
						+ (int) (windowSize / 2));
				min = min(quality, i - (int) (windowSize / 2), i
						+ (int) (windowSize / 2));
			}
			int actualWindowSize = 1
					+ Math.min(seqLength - 1, i + (int) (windowSize / 2))
					- Math.max(0, i - (int) (windowSize / 2));
			toReturn[i] = (double) sum(quality, i - (int) (windowSize / 2), i
					+ (int) (windowSize / 2))
					/ (actualWindowSize);
		}

		return toReturn;
	}

	public static int sum(byte[] values, int start, int stop) {
		start = Math.max(0, start);
		stop = Math.min(values.length, stop + 1);
		int sum = 0;
		for (int i = start; i < stop; i++) {
			sum += values[i];
		}
		return sum;
	}

	/**
	 * Returns the largest byte of the array in the subset start...stop.
	 * 
	 * @param values
	 *            data from which to pick largest.
	 * @param start
	 *            start position of subset to pick from.
	 * @param stop
	 *            stop position of the subset to pick from.
	 * @return
	 */
	public static byte max(byte[] values, int start, int stop) {
		start = Math.max(0, start);
		stop = Math.min(values.length, stop + 1);
		byte max = Byte.MIN_VALUE;
		for (int i = start; i < stop; i++) {
			if (values[i] > max) {
				max = values[i];
			}
		}
		return max;
	}

	/**
	 * Returns the smallest byte of the array in the subset start...stop.
	 * 
	 * @param values
	 *            data from which to pick smallest.
	 * @param start
	 *            start position of subset to pick from.
	 * @param stop
	 *            stop position of the subset to pick from.
	 * @return
	 */
	public static byte min(byte[] values, int start, int stop) {
		start = Math.max(0, start);
		stop = Math.min(values.length, stop + 1);
		byte min = Byte.MAX_VALUE;
		for (int i = start; i < stop; i++) {
			if (values[i] < min) {
				min = values[i];
			}
		}
		return min;
	}

}